##### ####################################################
#####                                               ######
#####             Commonly used functions             
#####                                               ######
##### ####################################################


block.bootstrap <- function(block){
  n.blocks <- length(unique(block))
  block.table <- table(block)
  resampled.blocks <- sample(unique(block),size=n.blocks,replace=TRUE)
  n.resamp <- sum(block.table[resampled.blocks])
  resample <- rep(NA,n.resamp)
  counter <- 0
  for (i in 1:n.blocks){
    tmp <- which(block == resampled.blocks[i])
    resample[counter + 1:length(tmp)] <- tmp
    counter <- counter + length(tmp)
  }
  return(resample)
}	


bootstrap.function <- function(model, cluster, data, reps=100, run.boot=T, weights = NULL){
  
  if(is.null(weights)) {
    data$model_weights <- 1
  }else{
    data$model_weights <- weights
  }
  
  ## Run the true model
  time <- system.time(true.model <- lm(model,data=data, weights = model_weights))
  
  ## Estimate time for bootstrapping
  print(paste("Estimated time:",round(as.numeric(time[3]*reps)),"seconds"))
  
  ## Initialise matrix to hold coefficients
  boot.coef <- matrix(NA, nrow=length(coef(true.model)), ncol=reps)
  rownames(boot.coef) <- names(coef(true.model))
  
  if(run.boot==T){
    
    ## Loop over reps
    for(r in 1:reps){
      cat(".")
      boot.data <- data[block.bootstrap(data[[cluster]])]
      boot.model <- lm(model, data=boot.data, weights = model_weights)
      boot.coef[names(coef(true.model))%in%names(coef(boot.model)),r] <- coef(boot.model) # This assigns coefficients to the correct elements in the boot matrix
    }
    
  }
  
  ##  Warning if coefficients are missing from some bootstrap samples
  if(any(is.na(boot.coef))) warning("Some coefficients are NA. Variances may be unreliable.")
  
  ## Save simulated coefficients
  
  true.model$simulations <- boot.coef
  
  ## Assign regular standard errors to the model object
  true.model$se <- summary(true.model)$coef[,2]
  
  ## Assign clustered standard errors to the model object
  true.model$clus.vcov <- vcovCluster(true.model,data[[cluster]])
  
  ## Assign clustered standard errors to the model object
  true.model$clustered.se <- coeftest(true.model,vcov=true.model$clus.vcov)[,2]
  
  ## Assign bootstrapped standard errors vector to the model object
  if(run.boot==T){
    true.model$boot.se <- apply(boot.coef,1,function(x) sd(x,na.rm=T))
  }else{
    true.model$boot.se	<-rep(NA, length(true.model$clustered.se))
  }
  
  ## Return the model object
  return(true.model)
}

mod_stargazer <- function(...){
  output <- capture.output(stargazer(...))
  # The first three lines are the ones we want to remove...
  output <- output[grep("begin\\{tabular\\}",output):(grep("textit\\{Note:\\}",output)-1)]
  # cat out the results - this is essentially just what stargazer does too
  cat(paste(output, collapse = "\n"), "\n")
}

# vcovCluster.r 
# function to compute var-cov matrix using clustered robust standard errors
# inputs:
# model = model object from call to lm or glm
# cluster = vector with cluster ID indicators
# output:
# cluster robust var-cov matrix
# to call this for a model directly use:
# coeftest(model,vcov = vcovCluster(model, cluster))
# formula is similar to Stata's , cluster command
#
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}


## ########
## Effect sizes
## ########

effect.func <- function(model, coefficient=2, baseline, se="boot.se"){
  
  effect <- coef(model)[coefficient]/baseline
  upper <- (coef(model)[coefficient]+1.96*model[[se]][coefficient])/baseline
  lower <- (coef(model)[coefficient]-1.96*model[[se]][coefficient])/baseline
  
  return(round(c(effect=as.numeric(effect), upper =as.numeric(upper), lower =as.numeric(lower))*100))
}

effect.func.gam <- function(model, coefficient=2, baseline, se="se"){
  
  effect <- coef(model)[coefficient]/baseline
  upper <- (coef(model)[coefficient]+1.96*summary(model)[[se]][coefficient])/baseline
  lower <- (coef(model)[coefficient]-1.96*summary(model)[[se]][coefficient])/baseline
  
  return(round(c(effect=as.numeric(effect), upper =as.numeric(upper), lower =as.numeric(lower))*100))
}


ci.calc <- function(model,coefficient=2){
  
  point <- coef(model)[coefficient]
  
  ## Regular standard errors
  upper <- point + 1.96 * model$se[coefficient]
  lower <- point - 1.96 * model$se[coefficient]
  
  ## Clustered standard errors
  upper.clustered <- point + 1.96 * model$clustered.se[coefficient]
  lower.clustered <- point - 1.96 * model$clustered.se[coefficient]	
  
  ## Bootstrapped standard errors
  upper.boot <- 2*point - quantile(model$simulations[coefficient,],c(0.025))
  lower.boot <- 2*point - quantile(model$simulations[coefficient,],c(0.975))	
  
  cis <- data.frame(rbind(
    c(coef=as.numeric(point), upper=as.numeric(upper), lower=as.numeric(lower)),
    c(coef=as.numeric(point), upper=as.numeric(upper.clustered), lower=as.numeric(lower.clustered)),
    c(coef=as.numeric(point), upper=as.numeric(upper.boot), lower=as.numeric(lower.boot))))
  
  cis$error <- c("Standard","Clustered","Block bootstrap")
  
  return(cis)
}


panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  txt <- paste(prefix, txt, sep="")
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}


collect_models_function <- function(data = rankings, dv = "eigen", bootstrap.reps=10, run.boot=T, formulas = list(), weights = NULL){
  
  if(length(formulas) == 0) stop("No formulas given. Needs to be a list of 'as.formula' objects.")
  
  model.out <- lapply(formulas, function(x){
    cat(".")
    this_mod <- as.formula(paste0(dv, x))
    mod <- bootstrap.function(this_mod, cluster = "debate_department", data = data, reps = bootstrap.reps, run.boot = run.boot, weights = weights)
    return(mod)
  })
  
  names(model.out) <- paste0("mod.",1:length(formulas))
  return(model.out)
}


leads <- function(x,nleads=1){
  return(	c(x[-c(1:nleads)],rep(0,nleads)))
}

lags <- function(x,nlags=1){
  return(c(rep(0,nlags),x[-c(((length(x)-nlags)+1):length(x))]))
}

# model <- mod1
# cluster <- "debate_department"
# data <- model.data
# model.type <- "ols"
cluster.function <- function(model, cluster, data, model.type = "ols", weights = NULL){
  
  # Apply weights
  if(is.null(weights)) {
    data$model_weights <- 1
  }else{
    data$model_weights <- weights
  }
  
  ## Run the true model
  if(model.type=="ols") true.model <- lm(model,data=data, weights = model_weights)
  if(model.type=="quasipoisson") true.model <- glm(model, data=data, family="quasipoisson")
  
  ## Save simulated coefficients
  
  ## Assign regular standard errors to the model object
  true.model$se <- summary(true.model)$coef[,2]
  
  ## Assign clustered standard errors to the model object
  
  if(is.null(true.model$na.action)){
    true.model$clus.vcov <- vcovCluster(true.model,data[[cluster]])
  }else{
    true.model$clus.vcov <- vcovCluster(true.model,data[[cluster]][-true.model$na.action]) 
  }
  true.model$clustered.se <- coeftest(true.model,vcov=true.model$clus.vcov)[,2]
  
  names(true.model$coefficients)[1] <- "Constant"
  names(true.model$se)[1] <- "Constant"
  names(true.model$clustered.se)[1] <- "Constant"
  
  names(true.model$coefficients)[grep("GenderF:minister_genderF",names(true.model$coefficients))] <- "Interaction"
  names(true.model$se)[grep("GenderF:minister_genderF",names(true.model$se))] <- "Interaction"
  names(true.model$clustered.se)[grep("GenderF:minister_genderF",names(true.model$clustered.se))] <- "Interaction"
  
  ## Return the model object
  return(true.model)
}

trim <- function(s) gsub("^[[:space:]]+|[[:space:]]+$","",s)
word.count <- function(str1) sapply(gregexpr("\\W+", str1), length) + 1
